#include "Constants.h"
#include <gsl/gsl_poly.h>
#include <math.h>
#include <vector>

using namespace std;

double BravaisIndex(double thD, double eta)
{
    double sinThD = sin(thD);
    return (sqrt(eta * eta - sinThD * sinThD) / cos(thD));
}

void Roots(double p, double etaPerpendicular, double phi, double * roots)
{
	double c = asin(1/etaPerpendicular);

	double A = -8 * p * c / (PI * PI * PI);
	double C = ( 6 * p * c / PI - 2);
	double D = p * PI - phi;

	if(abs(A) < EPS)
		roots[3] = gsl_poly_solve_quadratic(0, C, D, &roots[0], &roots[1]);
	else
		roots[3] = gsl_poly_solve_cubic(0, C/A, D/A, &roots[0], &roots[1], &roots[2]);
}

double FresnelParallel(double n2, double angle)
{
    double R = 1;
    double n1 = 1;
    double cos_gammaI = cos(angle);
    double a = ((n1 / n2) * sin(angle));
    double b = a * a;

    if (b > 1)
        return R;

    double cos_gammaT = sqrt(1 - b);

    R = (n2 * cos_gammaI - n1 * cos_gammaT) / 
		(n2 * cos_gammaI + n1 * cos_gammaT);

    return min(1, R * R);
}

double FresnelPerpendicular(double n2, double angle)
{
    double R = 1;
    double n1 = 1;
    double cos_gammaI = cos(angle);
    double a = ((n1 / n2) * sin(angle));
    double b = a * a;

    if (b > 1)
        return R;

    double cos_gammaT = sqrt(1 - b);

    R = (n1 * cos_gammaI - n2 * cos_gammaT) / 
		(n1 * cos_gammaI + n2 * cos_gammaT);

    return min(1, R * R);
}

double Fresnel(double etaPerpendicular, double etaParallel, 
			double angle)
{
    return 0.5f * (FresnelPerpendicular(etaPerpendicular, angle) + 
		FresnelParallel(etaParallel, angle));
}

double T(double absorption, double gammaT)
{
    //double l = 1 + cos(2 * gammaT);
	double l = cos(gammaT);
    return exp(-2 * absorption * l);
}

double A(double absorption, int p, double h, double refraction, 
			double etaPerpendicular, double etaParallel)
{
	double gammaI = asin(h);

	if(p == 0)
		return Fresnel(etaPerpendicular, etaParallel, gammaI);

	double gammaT = asin(h / etaPerpendicular);

	double fresnel = Fresnel(etaPerpendicular, etaParallel, gammaI);
	double invFresnel = Fresnel(1/etaPerpendicular, 1/etaParallel, gammaT);
	double t = T(absorption, gammaT);

	//return (1 - fresnel) * (1 - fresnel) * pow(invFresnel, p-1) * pow(t, p);
	return (1 - fresnel) * pow(invFresnel, p-1) * pow(t, p) * (1 - invFresnel);
}

double InverseFirstDerivate(double p, double etaPerpendicular, double h)
{
    double gammaI = asin(h);
    double c = asin( 1 / etaPerpendicular );
    double dGamma = (6 * p * c / PI - 2) - 
		3 * 8 * (p * c / (PI * PI * PI)) * gammaI * gammaI;

    return sqrt(1 - h * h) / dGamma;
}


double NP(int p, double phi, double thD)
{
	double refraction = eta;

	double etaPerpendicular = BravaisIndex(thD, refraction);
    double etaParallel = (refraction * refraction) / etaPerpendicular;	
	
	double roots[4];
	Roots(p, etaPerpendicular, phi, roots);

	double result = 0;

	for(int index = 0; index < roots[3]; ++index)
	{
		double gammaI = roots[index];

		double h = sin(gammaI);
		double finalAbsorption = A(absorption, p, h, refraction, 
			etaPerpendicular, etaParallel);
		double inverseDerivateAngle = InverseFirstDerivate(p, etaPerpendicular, h);

		result += finalAbsorption * 2 * abs(inverseDerivateAngle);
	}

	return min(1, result);
	
}
